#define _FILE_OFFSET_BITS 64
#include <stdio.h>
#include <common.c>

//since sub-abs-diff > 0.0, there must be one allele with direction != 0.
void get_direction(const unsigned int n_q[], const unsigned int t_q[], unsigned int n, int direction[]) {
    double n_sum = sum_array(n_q, n), t_sum = sum_array(t_q, n);
    unsigned int i;
    for(i = 0; i < n; i ++) {
        double diff = (t_q[i] / t_sum) - (n_q[i] / n_sum);
        if(diff > 0.0)
            direction[i] = 1;
        else if(diff < 0.0)
            direction[i] = -1;
        else
            direction[i] = 0;
    }
}
int one_allele_same_direction(const unsigned int n_q[], const unsigned int t_q[], unsigned int n, const int direction[]) {
    double n_sum = sum_array(n_q, n), t_sum = sum_array(t_q, n);
    unsigned int i;
    for(i = 0; i < n; i ++) {
        double diff = (t_q[i] / t_sum) - (n_q[i] / n_sum);
        if(direction[i] > 0 && diff > 0.0)
            return 1;
        else if(direction[i] < 0 && diff < 0.0)
            return 1;
    }
    return 0;
}

int main(int argc, char *argv[]) {
    FILE *n_f = fopen(argv[1], "r"), *t_f = fopen(argv[2], "r");

    unsigned int chr, pos;
    char n_info[ALLELE_INFO_LEN], t_info[ALLELE_INFO_LEN];
    while(to_same_site(n_f, t_f, &chr, &pos, n_info, t_info)) {
        unsigned char n_seq[SITE_ALLELE_NUM][ALLELE_LEN];
        unsigned int n_reads[SITE_ALLELE_NUM];
        unsigned int n_quality[SITE_ALLELE_NUM];
        unsigned int n_num = parse_allele_info(n_info, n_seq, n_reads, n_quality);

        unsigned char t_seq[SITE_ALLELE_NUM][ALLELE_LEN];
        unsigned int t_reads[SITE_ALLELE_NUM];
        unsigned int t_quality[SITE_ALLELE_NUM];
        unsigned int t_num = parse_allele_info(t_info, t_seq, t_reads, t_quality);

        unsigned int num = union_allele(n_seq, n_quality, n_num, t_seq, t_quality, t_num);
        double stat = sum_abs_diff_freq(n_quality, t_quality, num);
        if(stat <= 0.0)
            continue;
        int direction[SITE_ALLELE_NUM];
        get_direction(n_quality, t_quality, num, direction);//freq normal < tumor, is +1; -1; 0;

        #define TIMES 10000
        unsigned int i, more_times;
        for(i = 0, more_times = 0; i < TIMES; i ++) {
            random_value(n_quality, num);
            random_value(t_quality, num);
            if(one_allele_same_direction(n_quality, t_quality, num, direction) && sum_abs_diff_freq(n_quality, t_quality, num) >= stat)
                more_times ++;
        }
        fprintf(stdout, "%u\t%u\t%lf\t%g\t%s\t%s\n", chr, pos, stat, more_times / (double) TIMES, n_info, t_info);
    }
}//"./run normal.allele.qual1 tumor.allele.qual1 > somatic.txt"
